Convergence, Non- negativity and Stability of a New 
Milstein Scheme with Applications to Finance 

Desmond J. Higham* Xuerong Mao^ Lukasz Szpruch* 



Abstract 

We propose and analyse a new Milstein type scheme for simulating stochastic differential equations 
(SDEs) with highly nonlinear coefficients. Our work is motivated by the need to justify multi-level Monte 
Carlo simulations for mean-reverting financial models with polynomial growth in the diffusion term. We 
introduce a double implicit Milstein scheme and show that it possesses desirable properties. It converges 
strongly and preserves non-negativity for a rich family of financial models and can reproduce linear and 
nonlinear stability behaviour of the underlying SDE without severe restriction on the time step. Although 
the scheme is implicit, we point out examples of financial models where an explicit formula for the solution 
to the scheme can be found. 
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1 Introduction 

We study numerical approximation of the scalar stochastic differential equation (SDE) 

dx{t) = f(x(t))dt + g{x(t))dw(t). (f .1) 

Here x{t) € R for each t > 0, and, for simplicity, x(0) is taken to be constant. We assume that 
/ G C^R) and g G C 2 (R,R). Throughout we let (O, J", {J r t }t>o, R) be a complete probability space 
with a filtration {J~t}t>o satisfying the usual conditions, that is, right continuous and increasing while 
IFq contains all P-null sets, and we let w(t) be a Brownian motion defined on the probability space. 

Numerical approximations for equation are well studied in the case of globally Lipschitz continuous 
coefficients |17) . Super-linearly growing coefficients, however, raise many new questions. An important 
example is the Heston stochastic volatility 3/2-model [91 ITS] : 

dx(t) = x(t)(fi- ax(t))dt + Px(t) 3/2 dw(t), ^,a,/3>0. (1.2) 

This equation is also known as an inverse square-root process and was shown to be useful for modelling 
term structure dynamics Q] . We may also view (|1.2p as a stochastic extension to the logistic equation [JJ . 
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Recently [H] demonstrated that the standard Euler-Maruyama (EM) discretization scheme can diverge 
in strong and weak senses for SDEs with super- linear coefficients. However, in [26j[2Z] it is shown that an 
implicit Euler-type method strongly converges for nonlinearities similar to (|1.2I) . These positive results 
rely heavily on a one-sided Lipschitz condition satisfied by the drift coefficients of the SDE (jl.lj) , that 
is, for some constant K, 

(x-y)(f(x)-f(y))<K\x-y\ 2 , for x, y e R. (1.3) 

The solution of (ll.2j) is non-negative and condition (|1.3[) holds in the relevant region x, y > 0. However, 
in general, numerical approximations do not preserve non-negativity and hence convergence theorems 
developed in [H] cannot be applied in this situation. 

Preservation of positivity is a desirable modeling property, and, in many cases, non-negativity of the 
numerical approximation is needed in order for the scheme to be well defined. For example, evaluating the 
diffusion coefficient in (jl.2[) for a negative argument does not make sense. Many fixes have been proposed 
in the literature, but these can lead to substantial bias |20j . For more information about positivity 
preservation we refer the reader to [3J [TH (S3J HH] • It was shown in [TS] that a class of balanced methods 
can preserve positivity under an appropriate choice of the weight functions, but strong convergence was 
proved only under a global Lipschitz condition [2T]. Kahl et al. [TJ proved that the classical EM scheme 
does not preserve positivity for any scalar SDE. In the case of the drift implicit EM scheme positivity 
can be preserved, [28], if the drift coefficient has a very special form, for example as in the Ait-Sahalia 
interest rate model [2J. It was also shown in 14 that the drift implicit Milstein scheme applied to 

dx(t) = a(fi - x[t))dt + 0x{t) p dw(t) a, fi, /? > 0, and p e [0.5], (1.4) 

preserves positivity with no restriction on the time step if p — 0.5. For p E (0.5, 1] we need to restrict 
the step size to lie below f3~ 2 . 

Higher order approximation carries some pitfalls. It was demonstrated in [TO] that Milstein applied to 
a linear scalar SDE has worse stability properties than EM, even once we allow for implicitness in the 
drift coefficient. This is undesirable, particularly in the multi-level Monte Carlo (MLMC) setting [51 [5J, 
where we are required to use many simulations with large discretization time step. It is therefore natural 
to look for more advanced numerical techniques that automatically capture such a property. 

In order to address the issues mentioned above, we introduce a new (8, tr)-Milstein scheme for a general 
scalar SDE. Given any step size At, we define the partition V&t := {tk = kAt : k — 0, 1, 2, ...} of the half 
line [0, oo). Letting X tk denote the approximation to x(tk), with X to = x(0) and Awt k = w(tk+i) — w(tk), 
the (8, cr)-Milstein-scheme then has the following form 

X th+1 = X tk + 8f(X tk+1 )At + (1 - 8)f(X tk )At + g(X tk )Aw tk + ^(X^Awl 

(±_^ Llg{Xtk)At _ ° L i g{Xtk+i)At: (L5) 

where < 9, a < 1 are free parameters and L 1 = g-§^- We note that the (0, 0)-Milstein scheme reduces 
to classical Milstein [35]. We will sometimes refer to (1, 1)-Milstein as the double implicit scheme. The 
advantage of the extra degree of implicitness offered by er will become clear as we analyse the method. 
We note that the (8, cr)-Milstein scheme can be naturally derived from the Ito-Stratonovich expansion. 
Indeed, we can rewrite SDE (jl.ll) into its Stratonovich form 

dx{t) = f[x{t))dt + g(x(t)) o dw(t), 

where f(x) = f(x) — L 1 g(x). In the scalar case the drift-implicit Milstein scheme for the Stratonovich 
SDE is~given by (see [T7] p. 345) 

X tk+1 = X tk +l(X tk+1 )At + g(X tk )Aw tk + ii^JA^. 

Hence, we note that (8, cr)-Milstein may be obtained from the implicit Milstein scheme for a Stratonovich 
SDE. 
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In this work, we allow for a nonlinear drift coefficient and show that once p > 0.5 in (j 1 .4[) the step size 
restriction for non-negativity can be eliminated by the (9, cr)-Milstein method. We also present fairly 
general conditions for a family of Milstein schemes to preserve positivity. Due to that property the 
one-sided Lipschitz structure (| 1 .3[) will not be lost. Hence, the new scheme can be shown to converge 
strongly to the solution of the SDE (ll.2[) . Numerically we observe that the rate of strong convergence is 
1, which Giles [SJE] has shown to be the optimal rate from the MLMC perspective. 

The material is structured as follows. Section 2 contains proofs of the existence of positive solutions 
to (ll.5[) . In Section 3 we consider stability properties of the double implicit scheme. As motivation we 
demonstrate that for linear test SDEs we can significantly improve stability properties of the Milstein 
scheme. We then extend this result to a more general nonlinear setting. In Section 4 we develop the 
convergence results. We give conclusions in Section 5. 



2 Existence of a Solution for the Implicit Schemes 

We begin with conditions that guarantee the existence of a unique solution to (| 1 . 5|) . These will motivate 
the assumptions that we use to force positivity. 

Lemma 2.1. Let F be a function defined on R and consider the equation 

F(x) = b, (2.1) 
for a given i)£ R. If F is strictly monotone, i.e., 

(x - y)(F(x) - F(y)) > 0, (2.2) 
for all i,i/£R, x ^ y, then equation \2. 1\) has at most one solution. If F is continuous and coercive, 



i.e. 



lim ^£M = 0Oj (2.3) 

\x\^oo \X\ 

then for every b G M., equation &2.1)) has a solution Moreover, the inverse operator F~ 1 exists. 

A proof follows directly from Theorem 26. A in [29 . In order to prove that the (9, <r)-Milstein (|1.5[) 
scheme is well defined we impose two conditions. 

Assumption 2.2. Coefficients f and g in (jl.lj) are locally Lipschitz continuous and satisfy the following 
two conditions: 

One-sided Lipschitz condition. There exists a constant K > such that 

(x-y)(f(x)-f(y))<K\x-y\ 2 for all x,yeR. (2.4) 

Monotone condition. Operator L 1 acting on g satisfies 

(x-y)(L 1 g(x)-L 1 g(y))>0 for all x,yeR. (2.5) 

Remark 2.3. From Assumvtion \2.2\ and the Young inequality we may show that the drift coefficient f 
satisfies a one-sided Lipschitz-type condition 

xf(x) < K\x 2 \ + xf(0) < a + b \x\ for all x € R, 

where a = 0.5|/(0)| 2 and b — (2K + l)/2. Also from Assumvtion \2.2\ we can show that xL 1 g(x) is 
bounded below by a linear function 

xL 1 g{x) >xL 1 g(0) for all x,yeR. (2.6) 
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Lemma 2.4. Define, for any given At < (9(K + 1)) 1 , 

F{x) = x-9f(x)At + ^L 1 g(x), iel. (2.7) 

Then under Assumption for any b £ R there exists a unique x £ R such that Fix) — b and hence 
the method (|1.5p is well defined. 



Proof. In view of Lemma 1 2. II it is enough to show that the function F is continuous, coercive and strictly 
monotone. Clearly, F(x) is continuous on R. By Assumption 12.21 



(x - y)(F(x) - F(y)) >\x- y\ 2 - BKAt \x - y\ 2 = (1 - OK At) \x - y\ 2 > 0, 
for At < (9(K + Also by Assumption O and Remark [231 

xF(x) = x(x - 9f(x)At + ^L 1 g(x)At) (2.8) 
> \x\ 2 (1 - 2 -^±At) 6 -x |/(0)| 2 At + ^xL l g{Q)At, 
so F is coercive. □ 

From now on we assume that 

At < (9{K + I))- 1 (2.9) 



2.1 Existence of a Positive Solution for the (9, cr)-Milstein Scheme 

In this subsection we introduce assumptions on coefficients / and g of equation p. II) that allow us to 
prove the existence of a positive solution to (|1.5I) . 

Definition 2.5. Given x(0) > 0, if the solution of dTTTJ) satisfies ¥{{x(t) > : t > 0}) = 1 (¥({x(t) > 
: t > 0}) = 1), t/ien a stochastic one-step integration scheme computing approximations X tk ~ x(tk) 
preserves positivity (non-negativity) if 

F({X th+1 > 0\X tk > 0}) - 1 (P({A tfc+1 > 0\X tk > 0}) = 1). 

Let us note that to use the ideas from Lemma 12.41 to prove the existence of a positive solution to the 
implicit scheme we need to assume that a one-sided Lipschitz condition on / and monotone condition on 
L^g hold only on the positive domain. This significantly relaxes the conditions required for the existence 
and uniqueness of a solution to the implicit scheme (jl.5l) . 

Assumption 2.6. Coefficients f and g of the equation are locally Lipschitz continuous and satisfy 

the following two conditions: 

One-sided Lipschitz condition on R+ . There exists a constant K > such that 

(x-y)(f(x)- f(y))<K\x~y\ 2 for all x,y eR+. (2.10) 
Monotone condition on R + . Operator L 1 acting on g satisfies 

(.T-y)(L 1 .g( 2 ;)-L 1 ff (y))>0 for all x,y£R+. (2.11) 
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Many mean-reverting models with super- and sub-linear diffusion coefficients satisfy Assumption 12.61 for 
example, the mean-reverting SDEs 

dx(t) = (jj, - x(t) q )dt + x(t) p dw(t) for x G R, 

with fj,, q > and p > 0.5. 

In general, boundary behavior of one-dimensional SDEs can be fully characterized by the Feller test [IB] , 
Let us consider the interval (0,oo). We assume that / and g are locally Lipschitz continuous in (0,oo) 
and that g 2 (x) > 0, for x £ (0, oo). Let us also define the scale function 



p(x) = / exp 



9 2 {z) 



ds. 



where c G R. Since we analyse the behaviour of the above function at 0, we assume that c > x. By 
Proposition 5.22 in [TB] we have that if p(0+) = — oo then P [info<t<<x> x(t) = 0] = 1. Therefore, in order 
to show that the solution to is non-negative it is enough to show that p(0+) = — oo. 

Assumption 2.7. The coefficients f and g in satisfy the following conditions: 

f (0) > 0, g 2 {x) > for x e (0,oo), and g(Q) = for x = 0. 

To understand Assumption 12.71 better, we proceed with a heuristic argument. Lets assume that the 
solution of Ql.ljl attained at time t. Since the solution is Markovian we can consider the solution to 
this SDE with initial condition x(t) =0 that reads 

dx(t) = f(0)dt + g(0)dw(t). 
It is clear that we need to have g(0) = and /(0) > in order for x(t) to stay non-negative. 
Theorem 2.8. Let Assumvtions [^7E\ and \^7l\ hold. In addition we require that 

L 1 g(x) > for x > 0. 
Then there exists a unique positive solution to the {6,o~)-Milstein scheme (|1.5[) if 

x - + (! " 0)/(*) A * - i -^-L 1 g(x)At > -0/(O)At, x > 0. (2.12) 

2L L g{x) 2 

Similarly, a unique non-negative solution exists if 



2Li 5 (x) 



+ (1 — 6)f(x)At — - — - — -L g(x)At > —6f(0)At, x > 0. (2.13) 



Proof. In view of Lemma 12.11 and Definition 12.51 in order to prove the lemma we analyse the following 
equation 

F(X k+1 )=X k , Vfc. 

First we prove that P{A r / c+ i > | X k > 0} = 1. By Assumptions 12.61 operator F in (|2.7D is monotone on 
(0, oo) and we have 

,. xF(x) 

hm — j — j — = oo. 

By Assumption 12.71 we arrive at 

lim = -9f(0)At. 

x^0+ \x\ 

Hence operator F is coercive on (0, oo). Due to Lemma 12. 11 we may complete the proof by showing 
b(x) = x+(l-0)f(x)At + g(x)Aw tk+1+ ^L 1 g(x)Awl +1 - { ^-^L 1 g(x)At > -0f(O)At, for x > 0, 
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from which it follows that there exists a positive solution to F(X tk+1 ) = b(X tk ). First, for any given 
x > we find the minimum of the function 

H{y) = g{x)y + ^L 1 g{x)y 2 . 

Under the assumption L 1 g(x) > 0, for x > 0, this function possesses a global minimum 

m mH(y) = -^- 
y 2L L g(x) 

Hence 

b(x) > x + (1 - 9)}{x)At - ^—^L 1 g{x)At - -f^- > -0/(O)At, x > 0, 

2 lL L g[x) 

as required. For the non-negative case we have b(x) > —9f(0)At, x > 0. In that case we also need to 
check what happens if for some k we have the following event {Xk+i = | Xk > 0} (that corresponds to 
the case where b(x) = —9f(0)At ). Then by Assumption 12. 71 b(0) = (1 — 6)f(0)At and we require that 
b(0) > -9f(Q)At. That holds due to Assumption O □ 

For the fully implicit (1, 1)-Milstcin scheme we see from (|2.12[) that a condition guaranteeing non- 
negativity independently of At is 

> 0, x > 0. 



2Lig(x) 



2.2 Example: Heston Volatility Model 

Now we demonstrate that approximation of the 3/2-Heston volatility model (|1.2j) with the double implicit 
Milstein scheme preserves non-negativity. We point out that implicitness in the numerical approximation 
does not increase computational cost in this case, since we are able to find an explicit solution. This 
often will be the case in mathematical finance, where typical models have drift and diffusion coefficients 
of a polynomial type. 

The (1, 1)-Milstein scheme has the form 

X tk+1 =X tfe +/(X tfc+ JAi + 5 (X t JA Wtfc +iL\ 9 (X t JA<-iL 1 .g(X tfc+1 )At, (2.14) 

where now f(x) — [ix — ax 2 , g(x) = (3x 3 / 2 and L 1 g(x) = |/3 2 £ 2 . Clearly, the coefficients of equation (|1.2j) 
satisfy Assumptions [2~6l and [2~7l Hence, we may show that (|2 . 14[) has a unique non- negative solution by 
verifying condition (|2.13[) in Theorem 12.81 This reduces to x — x/3 > for x > and the result follows. 

An explicit formula for X tk+1 can be found by solving the relevant quadratic equation and choosing the 
positive solution, to give 

X tk+1 = Ma+^At))- 1 (J (1 - nAtf + 4(a + \p 2 )At{X tk + f3X?{ 2 Aw tk + 3/4^X^1)^(1-^ . 



3 Stability Analysis 

In this section we examine the global stability of the (a, 6')-Milstein scheme (jl.5D . The stability conditions 
we derive are related to mean-square stability, and we are interested in results that do not put severe 
restrictions on the time step. We begin with linear test equations where we can derive sharp results and 
represent stability regions graphically. 
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3.1 Linear Mean-Square Stability 

For the linear test SDE 

dx{t) = ax(t)dt + nx(t)dw(t), (3.1) 
the property of mean-square stability, 

lim E\x{t)\ 2 = 0, 

t— >oc 

is characterized by 

{2a + n 2 ) < 0. (3.2) 

For the #-Milstein scheme on (|3.1[) . 

Xt k+1 = X tk + 9aX tk+1 At + (1 - 6)aX tk At + ^X tk Aw tk+1 + \p?X tk [Aw 2 k+1 - At], 
the analogous property 

lim E|X t J 2 = 0, (3.3) 

k^oo 

was studied in [10 . In particular, the linear stability region 

Rms ■= {Ata, At// 2 e K : method mean-square stable on (j3.ip } (3-4) 

was shown to be significantly smaller than that for the corresponding Euler-based scheme. We now 
examine the new Milstein scheme (|1.5[) in this setting, which reduces to 

X tk+1 = X tk + eaX tk+1 At + (1 - 6)aX tk At + ^X tk Aw tk+1 (3.5) 

+ 1 r 2 X tk A w l +i ^l^ Xtk At-^X tk+1 At. 

Theorem 3.1. The (9, a) -Milstein scheme Q3.5[) is linearly mean-square stable, \3.S\) , if and only if 

(2a + p 2 ) + Ata 2 (l - 26) + ^-(2aa + /i 2 ) < 0. (3.6) 

Proof. We rewrite (|3. 51) as a recurrence of the form 

Xt k+1 =X tk (p + q£t k+1 +rg k+1 ) , 



where £ ~ N(0, 1) and 



1 + (1 - 6)aAt - t^l^At 
1 - 9aAt+£fi 2 At 

fiVAi 



1 -0aAt + §/i 2 Ai' 
i M 2 At 



Then 

\X I 2 



l-0oAt+f/x 2 At' 

\X tk | 2 (p 2 + q 2 & k+1 + r 2 tf k+1 + 2 Pq £ th+1 + 2 W e tk+1 + 2qr$ k+1 ] 
Taking conditional expectation of both sides lead us to 

n\Xt k+1 \ 2 \F th ] =\X tk \ (p 2 + q 2 + 3r 2 + 2pr) . 

Taking conditional expectation of both sides again we obtain 

E \X tk+1 1 2 = E \X tk | 2 (p 2 + q 2 + 3r 2 + 2pr) . (3.7) 

Therefore stability is characterized by (p + r) 2 + q 2 + 2r 2 < 1. This is equivalent to p.6[) . as required. □ 



Remark 3.2. Let us observe that for 8 = 0.5 and a = 1 we have recovered precisely the condition K3.2\\ 
for the underlying SDE, so the method perfectly reproduces stability for any step-size. More generally, 
for 8 > 0.5 and a — I we have the property that "problem stable implies method stable for all At", which 
is refered to as A-stability in the deterministic literature. 

Motivated by [10L 111] we will draw stability regions for (|3.5[) in the x-y plane, where x = aAt and 
y = /j, 2 At. In Figured] the stability region of the underyling SDE (I3.1[) is shaded light grey. The upper 
pictures in Figure [1] superimpose the stability region of the (8, 0)-Milstein scheme with 6 = 0,0.5,1, 
respectively, using darker shading. We see that even in the case of a linear scalar equation we are not 
able to reproduce the stability region of the underyling test equation (|3.1[) . However, by introducing 
additional implicitness we overcome this poor performance. The lower pictures in Figure [1] superimpose 
the stability region of the (8, cr)-Milstein scheme with (0, 1), (0.5, 1), (1, 1), respectively. As stated in 
Remark 13.21 we recover exactly the stability region of underlying test SDE (|3.ip for 8 = 0.5 and a = 1. 

9 = 0,o = 6 = 0.5,0 = 9= 1,0 = 




aAt aAt aAt 



Figure 1: Light shading: linear mean-square stability of the SDE. Darker shading: linear mean-square 
stability of Implicit Milstein (upper) and double-implicit Milstein (lower). 



3.2 Lyapunov Stability 

We begin this section by stating a result that combines Doob's Decomposition and Martingale Conver- 
gence Theorems. 

Theorem 3.3 (Lipster and Shiryaev |19j). Let Z = {Z„} n6 N be a nonnegative decomposable stochastic 
process with Doob-Meyer decomposition Z n = Zq + A„ — A„ + M n , where A 1 = {A^} ne N and A 1 = 
{A^}„gN are a.s. nondecreasing, predictable processes with Aq = Aq — 0, and M = {M n } nS N is local 
{J- n } n eN -martingale with Mq = 0. Then 

< U) : lim A l (n) < oo > C < u : lim A 2 (n) < oo > (1 \ lim Z n < oo > a.s. 
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In the authors proved a very general Stochastic LaSalle Theorem. Here we present a simplified 
version of their theorem, with fixed Lyapvmov function V(x) — \x\ . 

Theorem 3.4 (Shcn et al. [24 ). Let local Lipschitz conditions hold for f and g. Assume further that 
there exists a function z € C(R;K+) such that 



1 2 

xf{x) + - \g(x)\ < -z(x), 

for all x € R. For any xq G R, the solution (x(t)) t >o of (jl.ip then has the properties that 
limsup |a;(i)| 2 < oo a.s. and lim z(x(t)) = a.s. 



(3.8) 



Further if z(x) = if and only if x — 0, then 



lim x(t) = a.s. Vie 

t— >oo 



Now we present a counterpart of this Stochastic LaSalle Theorem for the new Milstein scheme. 

Theorem 3.5. Let Assumvtion \2. 2\ hold. Assume that for the (9, a) -Milstein Scheme (11.51) there exists 
a function z € C(K;R+) such that 



2xf(x) + \g(x)\ 2 + (l-26)\f(x)\ 2 At 
+ ^L 1 g(x)(2af(x) + L 1 g(x)) < -z(x) for all x £ 



(3.9) 



Then 



and 



lim sup \X tk | < oo 

k— >oo 



lim z(X tk ) = a.s. 



Further if z{x) = if and only if x — then 



lim X t , = a.s. 



Proof. We can rewrite F in (|2.7p as 

F(X tk+1 ) = F(X tk ) + f(X tk )At + g(X tk )Aw tk+1 + ^L 1 g(X tk )(Aw 2 tk+1 - At). 
Squaring both sides, we arrive at 

\F(X tk+1 )\ 2 = \F(X tk )\ 2 + \f(X tk )At\ 2 + \g(X tk )\ 2 At + \ \L x g{X th )f At 2 + 2F{X tk )f{X tk )At + m k+1 , 



where 



,,,,,, = \>, \ )\ 2 (Aw 2 k+1 - At) + \ (L^^JI 2 [(Aw 2 k+1 - At) 2 - 2At 2 } 

1 



g(X tk )Aw tk+1 + -i 1 . 9 (X t J(A< +i - At) 



+ 2F{X t 

+ 2f(X tk )At 
+ g(X tk )L 1 g(X tk )(Aw 2 k+i - At)Aw tk+1 



g{X tk )Aw tk+l + -i 1 . 9 (X t J(A< +i - At) 



(3.10) 



is a local martingale difference. From the definition of F we arrive at 

\F(X tk+1 )\ 2 = \F(X tk )\ 2 + 2X k f(X tk )At + \g(X tk )\ 2 At + ^gpQj [L^X^) + 2a/(X t J] At 2 

+ (1 - 29) \f{X tk )\ 2 At 2 + m k+1 . (3.11) 



Therefore 



where 



\F(X tk+1 )\ 2 = \F(X tk )\ 2 - A tk At + m k+1 , 



A k (x) = - ^2X tk f(X tk ) + \g(X tk )\ 2 + ^L 1 g(X tk )[L 1 g(X tk ) + 2af(X tk )] At 
+ (l-26)\f(X tk )\ 2 At). 
Hence, we have obtained a decomposition that allows us to apply Theorem 13.31 i.e., 

N N 

\F(X tN+1 )\ 2 = \F(X to )\ 2 -^A tfc At + ^m fc+1 . 

k=0 k=0 

Theorem [S3] gives lim*.-^ \F(X tk )\ 2 < oo. By condition fl3U) and (gUJ) 



1 x 2 



|F(a;)r = - 9f(x)At + -crL L g(x)At 

= \x\ 2 - 20xf{x)At - 9af{x)L 1 g{x)At 2 + 9 2 {f{x)) 2 At 2 + ^a 3 {L l g{x)f At 2 + axL 1 g{x)At 

> \x\ 2 + 9z(x)At + axL 1 g(x)At 

> \x\ 2 -a\x\ |L X ff(0)| At 

>(1 - 0.5At) |.t| 2 - 0.5cr 2 ^^(O)! 2 At. 

Hence lim sup j,.^ |X | 2 exists and is finite almost surely. Another implication of Theorem 13.31 is 

oo oo 

J2z(X tk )At<J2A tk At<oo a.s, 

k=0 k=0 

as required. □ 
In the case where (jl.5p is non-negative it is enough if condition p.9[) holds on the non-negative half line. 



Theorem 3.6. Let conditions required for existence of non-negative solution in Theorem ] 2. 8\ hold. As- 
sume that for the (9, o~)-Milstein Scheme (|1.5p there exists a function z 6 C(M;R+) such that 

2xf(x) + \g(x)\ 2 + (l-29)\f(x)\ 2 At 

+ ^L 1 g{x){2af(x) + L 1 ^)) < -z(x) for all x € R+. 



The 



limsup |X(£fc)| 2 < oo 

k— >oc 



lim z(X tk ) — a.s. 

fc— foo 



Further if z{x) = if and only if x = i/ien 



lim X tk — a.s. 

fc^oo 
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Proof. The proof is analogous to the proof of Theorem 13.51 □ 

Remark 3.7. Following on from Remark \3.Sl suppose that K3. 8}) holds, so the results of Theorem \3.4\ 
hold for the SDE. Then, to minimize restrictions on the stepsize in 13. 9\) . the choice 9 = 0.5 is clearly 
best, and the extra freedom allowed by the parameter a can be used to exploit dissipativity. For example, 
on the SDE 

dx(t) = -x{t) 3 dt + x(t) 2 dw(t), 

we have 

L 1 g(x)(2af(x) + L 1 g(x)) = L 1 g{x)(-2ax 3 + 2a; 3 ), 
so the choice a = 1 makes i3. 9]) independent of At and identical to h3.8\) . 



4 Convergence Result 



In this section we show that the numerical approximation (11.51) strongly converges to the solution of 
(|l.lj) under fairly general conditions. We will not establish the rate of convergence, but we perform 
numerical experiments that suggest a rate of 1. We note that the (1,1) scheme was considered in [T7] 
as an alternative to the more typical (1,0) version. In particular, those authors showed that when the 
coefficients /, g and L 1 g(x) in ([l.ip are globally Lipschitz, the (1, 1) case retains the usual first order of 
strong convergence. This result is easily extended to the general (9, a) case. 

Theorem 4.1. Let f, g and L 1 g(x) be globally Lipschitz. Then the (9,a)-Milstein scheme (|1.5[) strongly 
converges to the solution of the SDE (jl.ll) , that is 



E 



sup \x(t k ) 

0<t k <T 



X t k 



At p forp>2. 



Proof. A proof follows by extending the (1,1) case from Chapter 12 of Kloeden and Platen [T7]. □ 

Then using a localization procedure as in [5J 113] we can prove pathwise convergence without global 
Lipschitz Assumption. From [13] we know that scheme (jl.5l) almost surely converges to the solution of 
(TTTj) . that is: 

Theorem 4.2. Assume that the solution to SDE (|1.1[) has a strong solution. Then the {9 1 a)-Milstein 
scheme (jl.5l) converges to the solution of the SDE (II. ip in the pathwise sense, that is for 7 > there 
exists a random variable K — K(uj), uj € fi, such that 

sup \x(t k ) - X tk I < if(cj)At 1-7 , /orT, 7 >0. (4.1) 

0<t k <T 

Proof. A proof can be written in an analogous way to the proof of Theorem 1 in [13] . or Theorem 2.3 
in [8]. A key difference is that those authors used explicit schemes and defined continuous extensions 
to overcome some technical difficulties. In our setting, we need to define an J^-adapted continuous 
extension of the approximation (11.51) . Using notation of Lemma 2.4 we can write (ll.5[) in the form 

Xt k+1 = F-HX tk + (l-9)f(X tk )At + g(X tk )Aw tk + g(X tk )Aw 2 tk - ^^-L l g(X t k )At ) . 
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Then a suitable continuous extension of (|1.5[) for t € [t k ,t k+ i) could be defined by 

= ^ ^ + (1 - 0)f(X tk )(t - t fc ) + k )(w(t) ~ w(t k ) 
+ ^(^JM*) - w(t k f - ^l^g{X tk ){t - t k ^j . 

□ 

In order to show that we also have strong convergence we need to show that the solution to (| 1 . 5|) has 
bounded moments. We will prove boundedness of the moments under the following assumption. 

Assumption 4.3. Monotone-type condition. There exist constants a and b such that 
2xf(x) + \g(x)\ 2 + (l-26)\f(x)\ 2 At 

+ -^L 1 g(x)(2af(x) + L 1 g(x)) < a + b\x\ 2 forallxeR. (4.2) 

The following lemma establishes a useful relation between function F(x) defined in 12. 71 and its argument 
x. 



Lemma 4.4. Lets AssumytionsvLw and\4^ hold. Then there exist constants ci,C2 > such that 



\F{x)\ 2 > ci \x\ 2 - c 2 At forxe 



Proof. By Assumptions 12.21 and 14.31 we have 

\F(x)\ 2 > \x\ 2 - 8aAt - 9b \x\ 2 At + axL 1 g(x)At 

> \x\ 2 - 8aAt - 9b \x\ 2 At -a \x\ \L l g(0)\At 

> \x\ 2 - 6aAt - 6b \x\ 2 At - 0/2 \x\ 2 At - a 2 /(26) \L x g{ti)f At 

> (1 - (Ob + 6/2) At) \x\ 2 - 6aAt - a 2 /(26) (^(/(O)] 2 At, 



(4.3) 



and we take a = (1 - (6b + 6/2) At) and c 2 = -6a + a 2 / (26) \L 1 g(0)\ 2 . Due to ci > 0. □ 

Our analysis uses a localization procedure. We define the stopping time A m by 

X m =m£{k:\X tk \>m}. (4.4) 

We observe that when k £ [0, X m (u>)], \X tk _ 1 (uj)\ < m, but we might have that |X tfc (w)| > m, so the 
following lemma is not trivial. 



Lemma 4.5. Let Assumvtions [K^\ and \4-3\ hold. Then for p > 2 and sufficiently large integer m, there 
exists a constant K = K(p,m), such that 

E [\X tk \ p l [0 ,A m ](fc)] < K for any k > 0. 

Proof. By (|3.11[) and Assumption 14.31 we obtain 

\F(X tk )\ 2 < \F(X(t k ^)\ 2 +aAt + b \Xt k ^\ 2 At + Am k , 
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where Arrik+i is defined by (|3.10[) . Using the basic inequality (a\ + 0,2 + a 3 + a 4 ) p / 2 < 4 p / 2 1 (a^ + + 
a 3 + a 4:)y where > 0, we obtain 

,p/2\ 



\F(X tk )\ P <4 p - X ({FiX^ \ p + (aAtf/z + b \X tk _, \ p At + \Am k \* 



(4.5) 



As a consequence 



E[\F(X tk )\ p l [0 , Xm] (k)] <4P-^ (E[|F(X tfe _ 1 | p l [0 , Am] (A:)] + (aAt)^ 



■ bm p At + E 



\Am k \ p / 2 l [0tXm] (k) 



In order to bound E 



\Am,k\ p ^ 2 lm a i(fc) we need to consider all the terms of Am^ separately. By the 



Cauchy-Schwarz inequality 



E 



< 



AW. - At 



(E[|g(I tt J| 2p l [w W]) 1/2 (E 



p/2 



1 [0,A m ](fc) 



1/2. 



Aiu, - At 



I 2p 

Since there exists a positive constant C(p), such that E Aw^ J < C(p), there exists a constant 
C{m,p) such that 



P \l/2 



E 



Aw? - At 



l[o,A m ](*) < C(m,p). 



In the same way we can bound all the other terms of Am k . Hence 

E[|F(X t J| p l [0 , Am] (fc)] <C(m,p) 
Due to Lemma 14^41 the proof is complete. 



□ 



In addition to Assumption 14.31 we require the following very mild restriction on the coefficients of the 
SDE. 

Assumption 4.6. The coefficients of equation M.l\) satisfy a polynomial growth condition. That is, 
there exists a pair of constants h > 1 and H > such that 



\f(x)\V\g{x)\<H(l + \x\ h ), Vs. 



(4.6) 



Now we formulate the key theorem that allows us to prove a strong convergence result. 



Theorem 4.7. Let Assunwtions Wl^. |^.ff| and \4-6] hold. Then there exists a constant K — K(T) such 
that the (9,a)-Milstein scheme (|1.5p satisfies 

sup E\X tk \ 2 <K. 
a<t k <T 



Proof. By (|3.11D and Assumption 14. 31 we arrive at 



\F(X tk+ X < \F(X tk )\ 2 +aAt + b \X tk \ 2 At + Am k+1 , 



(4.7) 
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where Am,k+i is defined by (|3.10[) . Let N be any non-negative integer such that NAt < T. Summing 
both sides of inequality (|4.7[) from k — to N A A m , we get 

NA\ m NAXrn 

\F(X tNAXm+1 )\ 2 <\F(X t0 )\ 2 + aT + b \Xt k \ 2 At+ £ Am k+1 



fc=0 
N 



k=0 
N 



< \F(X t0 )\ 2 + aT + bJ2 \Xt kAXm \ Ai + Am fe+ il [0)Am] (fe). 



fe=0 



k=0 



Due to Lemma |4~5"1 T^^Lq Anih+i l[o ; A m ] (fc) i s a martingale. Hence 



E|^(X tNAAm+1 )| < |F(X t0 )|+aT + 6E 
Due to Lemma l4~4l we have 



,fc=l 



At 



E |F(X tjVAAm+1 )| 2 < |F(X t0 )| + (a + c 2 q 1 ) T + 6^1 
By the discrete Gronwall Lemma 



/V 



.fe=0 



At 



E |F(X tNAAm+1 )| 2 < [|F(X t0 )| + (a + c 2 q 1 ) T] exp (bc^ 1 T) 



(4.8) 



where we used the fact that NAt < T. Thus, letting m — > oo in (|4.8p and applying Fatou's lemma, we 
obtain 

E\F(X tN+1 )\ 2 < [\F(X t0 )\ + (a + c 2 c^)T] exp {bc^T). 
The final bound follows from Lemma |4~41 □ 



We are ready to prove a strong convergence result. 

Theorem 4.8. Let Assumvtions \2.2\ \4-S\ and \4-6\ hold. Then the (9,a)-Milstein scheme (11.51) strongly 
converges to the solution of the SDE (jl.ip , that is 



lim E \x(t k ) - X tk \ p = for < p < 2. 



(4.9) 



Proof. By (|4.1I) the (6, cr)-Milstein approximnation (|1.5I) X tfe converges to x(tfe) in probability (Theorem 
2.2 in [25 ). Theorem 14.71 implies that the sequence {|^t fc | 2 ~ £ }t fc is uniformly integrable (Lemma 2.3 in 
[25])- Therefore by the Vitali convergence theorem (Theorem 2.4 in 25 ) the statement of the theorem 
holds. □ 



In case where we can guarantee non-negativity of approximation, conditions required to prove Theorem 
IB] can be significantly relaxed. 



(4-10) 



Assumption 4.9. Monotone-type condition on K + . There exist constants a and b such that 
2xf(x) + \g(x)\ 2 + (l-20)\f(x)\ 2 At 

+ -^-L 1 g(x)(2af(x) + L 1 g{x)) <a + b \x\ 2 for all x € R+. 

Theorem 4.10. Let conditions required for existence of non-negative solution in Theorem \2.8\ hold. 
Then under Assumptions ^. 6] and \4-9\ the (9,o~)-Milstein scheme (11.51) strongly converges to the solution 
of the SDE CLU), that is 

A lim o E|a;(t fe )-X t J p = for < p < 2. (4.11) 
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Proof. The theorem can be proved in an analogous way to Theorem 14.81 



□ 



It is clear that 3/2-model (|1.2p doest not satisfy Assumption 14.31 but satisfies Assumption 14.91 as long as 
a > /? 2 /4. These condition seems not be restrictive as pointed out in [7]. 



4.1 Numerical Experiment 



In order to estimate the rate of convergence we proceed with numerical experiments for (|1.2p . We focus 
on the strong endpoint error, e^ ong = E \x(T) - X T \, with T = 1. We used fi = 0.1, a = 0.2, (3 = 
and x(0) = 0.5. We plot e^ ons against At on a log-log scale. Error bars representing 95% confidence 
intervals are shown by circles, and a reference line of slope 1 is also given. Although we do not know 




Figure 2: Strong error of double-implicit Milstein scheme applied to Heston 3/2 Stochastic volatility 
model. 

the explicit form of the solution, Theorem 14.81 guarantees that the (l.l)-Milstein scheme (|2.14l) strongly 
converges to the true solution. We therefore take the (l,l)-Milstein scheme with At = 2~ 14 as a reference 
solution. We compare this with the (l,l)-Milstein scheme evaluated with (2At, 2 3 Ai, 2 5 At, 2 7 At) in order 
to estimate the rate of convergence. Since we are using a Monte Carlo method, the sampling error decays 
like 1/vMj where M = 10000 is the number of sample paths. From Figure [2] we see that there appears 
to exist a positive constant C such that 

e At° ng — CAt, for sufficiently small At. 

A least squares fit for equality produced the value 1.1304 for the rate with residual of 0.2468. Hence, 
our results are consistent with strong order of convergence equal to one. 
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5 Conclusions 



Our aim was to introduce a new discretization scheme that can be shown to work well on highly nonlinear 
SDEs arising in mathematical finance and to possess excellent linear and nonlinear stability properties. 
There are several interesting areas for follow-up work; most notably (a) establishing a strong order of 
convergence for this method in a nonlinear setting, and (b) developing a theory of positivity preservation 
in the case of SDE systems and their numerical simulation. 
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